Drawing

Real-time $R_t$ of COVID-19 pandemic in Texas (Until Sep 22)¶

This report estimates the county level real-time effective production number ($R_t$) of the COVID-19 pandemic in Texas, based on data from the Texas Health and Human Services. This report is updated daily with the most recent version available at the ICTR website.

This report is created and maintained by data scientists from the Institute for Clinical & Translational Research, Baylor College of Medicine, with an algorithm that is based on Kevin Systrom's work for the US data. Please refer to the original report for details.

county_name = 'Harris'

def prepare_cases(cases, cutoff=24, smooth_days=7, smooth_std=2):
    new_cases = cases.diff()
    
    smoothed = new_cases.rolling(smooth_days,
        win_type='gaussian',
        min_periods=1,
        center=True).mean(std=smooth_std).round()
    
    idx_start = np.searchsorted(smoothed, cutoff)
    
    smoothed = smoothed.iloc[idx_start:]
    original = new_cases.loc[smoothed.index]
    
    return original, smoothed

def plot_cases(county_name, ax, cutoff=1):
    cases = counties.xs(county_name).rename(f"{county_name} cases")
    original, smoothed = prepare_cases(cases, cutoff=cutoff)
    
    index = original.index.get_level_values('date')
    values = original.values
    ax.plot(index, values, 
               c='k',
               linestyle=':',
               alpha=.5,
               label='Actual')
    ax.set_title(f"{county_name} New Cases per Day")
    index = smoothed.index.get_level_values('date')
    values = smoothed.values
    ax.scatter(index, values, label='Smoothed', lw=0.3, s=20)
    ax.get_figure().set_facecolor('w')
    #ax.set_xticklabels(, rotation=45, ha='right')
    
    ax.xaxis.set_major_locator(mdates.MonthLocator())
    ax.xaxis.set_major_formatter(mdates.DateFormatter('%b'))
    ax.xaxis.set_minor_locator(mdates.DayLocator())    
    ax.set_xlim(pd.Timestamp('2020-03-01'), original.index.get_level_values('date')[-1]+pd.Timedelta(days=1))

#county_name = 'Brewster'
cases = counties.xs(county_name).rename(f"{county_name} cases")
original, smoothed = prepare_cases(cases)
#fig, ax = plt.subplots(figsize=(600/72,400/72))
#plot_cases(county_name, ax)

def get_posteriors(sr, sigma=0.15):

# (1) Calculate Lambda
    lam = sr[:-1].values * np.exp(GAMMA * (r_t_range[:, None] - 1))

# (2) Calculate each day's likelihood
    likelihoods = pd.DataFrame(
        data = sps.poisson.pmf(sr[1:].values, lam),
        index = r_t_range,
        columns = sr.index[1:])
    
    # (3) Create the Gaussian Matrix
    process_matrix = sps.norm(loc=r_t_range,
                              scale=sigma
                             ).pdf(r_t_range[:, None])

# (3a) Normalize all rows to sum to 1
    process_matrix /= process_matrix.sum(axis=0)
    
    # (4) Calculate the initial prior
    #prior0 = sps.gamma(a=4).pdf(r_t_range)
    prior0 = np.ones_like(r_t_range)/len(r_t_range)
    prior0 /= prior0.sum()

# Create a DataFrame that will hold our posteriors for each day
    # Insert our prior as the first posterior.
    posteriors = pd.DataFrame(
        index=r_t_range,
        columns=sr.index,
        data={sr.index[0]: prior0}
    )
    
    # We said we'd keep track of the sum of the log of the probability
    # of the data for maximum likelihood calculation.
    log_likelihood = 0.0

# (5) Iteratively apply Bayes' rule
    for previous_day, current_day in zip(sr.index[:-1], sr.index[1:]):

#(5a) Calculate the new prior
        current_prior = process_matrix @ posteriors[previous_day]
        
        #(5b) Calculate the numerator of Bayes' Rule: P(k|R_t)P(R_t)
        numerator = likelihoods[current_day] * current_prior
        
        #(5c) Calcluate the denominator of Bayes' Rule P(k)
        denominator = np.sum(numerator)
        
        # Execute full Bayes' Rule
        posteriors[current_day] = numerator/denominator
        
        # Add to the running sum of log likelihoods
        log_likelihood += np.log(denominator)
    
    return posteriors, log_likelihood

# Note that we're fixing sigma to a value just for the example
posteriors, log_likelihood = get_posteriors(smoothed, sigma=.25)

Plotting in the Time Domain with Credible Intervals¶

Since our results include uncertainty, we'd like to be able to view the most likely value of $R_t$ along with its highest-density interval.

def plot_rt(result, ax, county_name):
    
    
    # Colors
    ABOVE = [1,0,0]
    MIDDLE = [1,1,1]
    BELOW = [0,0,0]
    cmap = ListedColormap(np.r_[
        np.linspace(BELOW,MIDDLE,25),
        np.linspace(MIDDLE,ABOVE,25)
    ])
    color_mapped = lambda y: np.clip(y, .5, 1.5)-.5
    
    index = result['ML'].index.get_level_values('date')
    values = result['ML'].values
    ax.set_title(f'{county_name} Real Time Rt ({values[-1]:.1f} on {result.index.get_level_values("date")[-1].strftime("%m/%d")})')
    
    # Plot dots and line
    ax.plot(index, values, c='k', zorder=1, alpha=.25)
    ax.scatter(index,
               values,
               s=40,
               lw=.5,
               c=cmap(color_mapped(values)),
               edgecolors='k', zorder=2)
    
    # Aesthetically, extrapolate credible interval by 1 day either side
    lowfn = interp1d(date2num(index),
                     result['Low_90'].values,
                     bounds_error=False,
                     fill_value='extrapolate')
    
    highfn = interp1d(date2num(index),
                      result['High_90'].values,
                      bounds_error=False,
                      fill_value='extrapolate')
    
    extended = pd.date_range(start=pd.Timestamp('2020-03-01'),
                             end=index[-1]+pd.Timedelta(days=1))
    
    ax.fill_between(extended,
                    lowfn(date2num(extended)),
                    highfn(date2num(extended)),
                    color='k',
                    alpha=.1,
                    lw=0,
                    zorder=3)

ax.axhline(1.0, c='k', lw=1, label='$R_t=1.0$', alpha=.25);
    
    # Formatting
    ax.xaxis.set_major_locator(mdates.MonthLocator())
    ax.xaxis.set_major_formatter(mdates.DateFormatter('%b'))
    ax.xaxis.set_minor_locator(mdates.DayLocator())
    
    ax.yaxis.set_major_locator(ticker.MultipleLocator(1))
    ax.yaxis.set_major_formatter(ticker.StrMethodFormatter("{x:.1f}"))
    ax.yaxis.tick_right()
    ax.spines['left'].set_visible(False)
    ax.spines['bottom'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.margins(0)
    ax.grid(which='major', axis='y', c='k', alpha=.1, zorder=-2)
    ax.margins(0)
    ax.set_ylim(0.0, 5.0)
    ax.set_xlim(pd.Timestamp('2020-03-01'), result.index.get_level_values('date')[-1]+pd.Timedelta(days=1))
    fig.set_facecolor('w')

#fig, ax = plt.subplots(figsize=(600/72,400/72))

#plot_rt(result, ax, county_name)
#ax.set_title(f'Real-time $R_t$ for {county_name} (till {result.index.get_level_values("date")[-1].strftime("%m/%d")})')
#ax.xaxis.set_major_locator(mdates.WeekdayLocator())
#ax.xaxis.set_major_formatter(mdates.DateFormatter('%b %d'))
#plt.xticks(rotation=45)

sigmas = np.linspace(1/20, 1, 20)

targets = ~counties.index.get_level_values('county').isin(FILTERED_REGION_CODES)
counties_to_process = counties.loc[targets]

results = {}

for county_name, cases in counties_to_process.groupby(level='county'):
    
    print(county_name)
#(county_name == "Coryell")|(county_name == "Jones")|(county_name == "Moore")|(county_name == "Panola")|(county_name == "Taylor"):
    if (county_name == "Bowie")|(county_name == "Coryell")|(county_name == "Harrison")|(county_name == "Houston")|(county_name == "Jones")|(county_name == "Panola")|(county_name == "Pecos")|(county_name == "Taylor"):
        continue
        
    new, smoothed = prepare_cases(cases, cutoff=1)
      
    max_cases = smoothed.max()
    if np.isnan(max_cases) or max_cases <= 2:
        continue
    
    if len(smoothed) == 0:
        new, smoothed = prepare_cases(cases, cutoff=1)
    if len(smoothed) == 0:
        continue
    result = {}
    
    # Holds all posteriors with every given value of sigma
    result['posteriors'] = []
    
    # Holds the log likelihood across all k for each value of sigma
    result['log_likelihoods'] = []
    
    for sigma in sigmas:
        posteriors, log_likelihood = get_posteriors(smoothed, sigma=sigma)
        result['posteriors'].append(posteriors)
        result['log_likelihoods'].append(log_likelihood)
    
    # Store all results keyed off of county name
    results[county_name] = result
    clear_output(wait=True)

print('Done.')

Daily new cases and real time $R_t$ for Texas Medical Center¶

def plot_tmc_cases(cases, ax, cutoff=15, smooth_days=7, smooth_std=2):
    original, smoothed = prepare_cases(cases, cutoff=cutoff, smooth_days=smooth_days,
                                           smooth_std=smooth_std)
    
    index = original.index.get_level_values('Collected Date')
    values = original.values
    ax.plot(index, values, 
               c='k',
               linestyle=':',
               alpha=.5,
               label='Actual')
    ax.set_title(f"TMC New Cases per Day")
    index = smoothed.index.get_level_values('Collected Date')
    values = smoothed.values
    ax.scatter(index, values, label='Smoothed', lw=0.3, s=20)
    ax.get_figure().set_facecolor('w')
    #ax.set_xticklabels(, rotation=45, ha='right')
    
    ax.xaxis.set_major_locator(mdates.MonthLocator())
    ax.xaxis.set_major_formatter(mdates.DateFormatter('%b'))
    ax.xaxis.set_minor_locator(mdates.DayLocator())    
    ax.set_xlim(pd.Timestamp('2020-03-01'), original.index.get_level_values('Collected Date')[-1]+pd.Timedelta(days=1))
    
    
def plot_tmc_rt(result, ax):
    # Colors
    ABOVE = [1,0,0]
    MIDDLE = [1,1,1]
    BELOW = [0,0,0]
    cmap = ListedColormap(np.r_[
        np.linspace(BELOW,MIDDLE,25),
        np.linspace(MIDDLE,ABOVE,25)
    ])
    color_mapped = lambda y: np.clip(y, .5, 1.5)-.5
    
    index = result['ML'].index.get_level_values('Collected Date')
    values = result['ML'].values
    ax.set_title(f'TMC Real Time Rt ({values[-1]:.1f} on {result.index.get_level_values("Collected Date")[-1].strftime("%m/%d")})')
    
    # Plot dots and line
    ax.plot(index, values, c='k', zorder=1, alpha=.25)
    ax.scatter(index,
               values,
               s=40,
               lw=.5,
               c=cmap(color_mapped(values)),
               edgecolors='k', zorder=2)
    
    # Aesthetically, extrapolate credible interval by 1 day either side
    lowfn = interp1d(date2num(index),
                     result['Low_90'].values,
                     bounds_error=False,
                     fill_value='extrapolate')
    
    highfn = interp1d(date2num(index),
                      result['High_90'].values,
                      bounds_error=False,
                      fill_value='extrapolate')
    
    extended = pd.date_range(start=pd.Timestamp('2020-03-01'),
                             end=index[-1]+pd.Timedelta(days=1))
    
    ax.fill_between(extended,
                    lowfn(date2num(extended)),
                    highfn(date2num(extended)),
                    color='k',
                    alpha=.1,
                    lw=0,
                    zorder=3)

ax.axhline(1.0, c='k', lw=1, label='$R_t=1.0$', alpha=.25);
    
    # Formatting
    ax.xaxis.set_major_locator(mdates.MonthLocator())
    ax.xaxis.set_major_formatter(mdates.DateFormatter('%b'))
    ax.xaxis.set_minor_locator(mdates.DayLocator())
    
    ax.yaxis.set_major_locator(ticker.MultipleLocator(1))
    ax.yaxis.set_major_formatter(ticker.StrMethodFormatter("{x:.1f}"))
    ax.yaxis.tick_right()
    ax.spines['left'].set_visible(False)
    ax.spines['bottom'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.margins(0)
    ax.grid(which='major', axis='y', c='k', alpha=.1, zorder=-2)
    ax.margins(0)
    ax.set_ylim(0.0, 5.0)
    ax.set_xlim(pd.Timestamp('2020-03-01'), 
                result.index.get_level_values('Collected Date')[-1] + pd.Timedelta(days=1))
    fig.set_facecolor('w')

Here are estimated $R_t$ in the past 15 days:

	ML	Low_90	High_90
Collected Date
2020-09-07	0.77	0.40	1.12
2020-09-08	0.82	0.43	1.16
2020-09-09	0.86	0.49	1.22
2020-09-10	0.89	0.49	1.23
2020-09-11	0.89	0.49	1.23
2020-09-12	0.87	0.47	1.21
2020-09-13	0.86	0.47	1.21
2020-09-14	0.86	0.45	1.20
2020-09-15	0.87	0.47	1.22
2020-09-16	0.87	0.48	1.23
2020-09-17	0.83	0.42	1.18
2020-09-18	0.77	0.37	1.13
2020-09-19	0.69	0.29	1.06
2020-09-20	0.60	0.20	0.98
2020-09-21	0.52	0.11	0.89

Daily new cases and real time $R_t$ for Harris County¶

NOTE: New cases reported on September 16 (2768) includes cases in the past month that were not reported in time. The estimated real-time $R_t$ for September 16 is therefore inaccurate.

Here are estimated $R_t$ in the past 15 days:

		ML	Low_90	High_90	Low_50	High_50
county	date
Harris	2020-09-07	0.65	0.44	0.80	0.54	0.69
	2020-09-08	0.59	0.40	0.76	0.50	0.65
	2020-09-09	0.68	0.48	0.85	0.60	0.75
	2020-09-10	0.80	0.61	0.98	0.69	0.85
	2020-09-11	0.99	0.79	1.16	0.87	1.03
	2020-09-12	1.21	1.00	1.37	1.12	1.27
	2020-09-13	1.48	1.28	1.64	1.38	1.53
	2020-09-14	1.56	1.38	1.72	1.48	1.62
	2020-09-15	1.58	1.39	1.72	1.48	1.62
	2020-09-16	1.43	1.25	1.57	1.35	1.48
	2020-09-17	1.24	1.05	1.37	1.15	1.28
	2020-09-18	0.96	0.78	1.10	0.88	1.01
	2020-09-19	0.77	0.60	0.92	0.66	0.80
	2020-09-20	0.59	0.42	0.75	0.49	0.63
	2020-09-21	0.65	0.47	0.81	0.57	0.71

Daily new cases and real time $R_t$ for all Texas counties¶

Most recent $R_t$ by county¶

filtered = final_results.index.get_level_values(0).isin(FILTERED_REGION_CODES)
mr = final_results.loc[~filtered].groupby(level=0)[['ML', 'High_90', 'Low_90']].last()

def plot_standings(mr, figsize=None, title='Most Recent $R_t$ by county'):
    if not figsize:
        figsize = ((15.9/50)*len(mr)+.1,2.5)
        
    fig, ax = plt.subplots(figsize=figsize)

ax.set_title(title)
    err = mr[['Low_90', 'High_90']].sub(mr['ML'], axis=0).abs()
    bars = ax.bar(mr.index,
                  mr['ML'],
                  width=.825,
                  color=FULL_COLOR,
                  ecolor=ERROR_BAR_COLOR,
                  capsize=2,
                  error_kw={'alpha':.5, 'lw':1},
                  yerr=err.values.T)

#     for bar, county_name in zip(bars, mr.index):
#         if county_name in no_lockdown:
#             bar.set_color(NONE_COLOR)
#         if county_name in partial_lockdown:
#             bar.set_color(PARTIAL_COLOR)

labels = mr.index.to_series().replace({'District of Columbia':'DC'})
    ax.set_xticklabels(labels, rotation=90, fontsize=11)
    ax.margins(0)
    ax.set_ylim(0, 5.)
    ax.axhline(1.0, linestyle=':', color='k', lw=1)

#     leg = ax.legend(handles=[
#                         Patch(label='Full', color=FULL_COLOR),
#                         Patch(label='Partial', color=PARTIAL_COLOR),
#                         Patch(label='None', color=NONE_COLOR)
#                     ],
#                     title='Lockdown',
#                     ncol=3,
#                     loc='upper left',
#                     columnspacing=.75,
#                     handletextpad=.5,
#                     handlelength=1)

#     leg._legend_box.align = "left"
    fig.set_facecolor('w')
    return fig, ax

mr.sort_values('ML', inplace=True)
plot_standings(mr);